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ABSTRACT 

This paper describes SEDfit, the earliest — but continually upgraded — software package for 
spectral energy distribution fitting (SED fitting) of high-redshift photometric data, and the only one 
to properly treat non-detections. The principles of maximum-likelihood SED fitting are described, 
including formulae used for fitting both detected and un-detected (upper limits) photometric data. 
The internal mechanics of the SEDfit package are presented and several illustrative examples of its 
use are given. The paper concludes with a discussion of several issues and caveats applicable to 
SED-fitting in general. 

Subject headings: methods: data analysis — techniques: photometric — galaxies: fundamental pa- 
rameters — galaxies: high redshift — galaxies: stellar content 
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1. INTRODUCTION 

The study of how galaxies form and evolve is one of the 
most active fields of astrophysical research today. Key 
in this work is the ability to estimate the stellar masses, 
star formation rates (SFRs), and dust content of faint 
galaxies in the distant Universe. 

One of the most widely-used approaches to determin- 
ing these quantities has been through the quantitative 
comparison of observed multi- wavelength spectral energy 
distributions (SEDs) with theoretical or empirical mod- 
els. This approach was first developed for high- 2 galax- 
ies by Sawicki & Yee (1998) and has subsequently been 
adopted and elaborated on by numerous other studies 
(e.g., Ouchi, Yamada, Kawai, & Ohta 1999; Papovich 
et al. 2001; Hall et al. 2001; Shapley et al. 2001; and, 
subsequently, many others). It is fair to say that SED fit- 
ting has now become a standard way to estimate stellar 
masses and other properties of high-redshift galaxies. 

A related task is the estimation of redshift from broad- 
band photometry, a technique known as photometric 
rcdshifts (photo-z's). When photo- z estimation is done 
by means of the template fitting approach (Gwyn & 
Hartwick 1996; Sawicki, Lin, & Yee 1997), it can be re- 
garded as a subset of the information that can be rou- 
tinely extracted from full SED-fitting. Indeed, photo- z 
estimation by template fitting has led to the subsequent 
development of SED-fitting as a tool for high-redshift 
galaxy studies. 

This paper gives a description of the SEDfit software 
tool, which is the first — though continually-evolving — 
among many software packages now used in SED-fitting 
work; it is also the only such package at present to prop- 
erly treat upper limits — i.e., observed but undetected 
photometric data. To date, SEDfit has been used to es- 
timate or simulate photometric redshifts (Sawicki, Lin, 
& Yee, 1997; Hogg et al. 1998; Sawicki 2002; Sawicki 
& Mallen-Ornelas 2003; Sorba & Sawicki 2010; Sorba & 
Sawicki 2011), to determine the physical properties of 
distant galaxies such as their stellar mass, starburst age, 
and reddening (Sawicki & Yee 1998; Sawicki 2001; Yabe 
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et al. 2009; Yuma et al. 2010, 2011; Hashimoto et al. 
2010; Palmer 2010, 2012), or to do both at the same time 
(Thompson et al. 1999; Hall et al., 2001; Sawicki 2012). 
It has also been used for simpler, but nevertheless impor- 
tant, tasks such as calculating colors and k-corrections of 
distant galaxies (Sawicki et al. 2005; Sawicki & Thomp- 
son 2006a, 2006b, Arcila-Osejo, Sawicki, & Sato 2012). 
It continues to be used in these and similar roles. 

The paper is organized as follows. Section [2] presents 
an overview of SEDfit, while subsequent sections delve 
into some of the key details: $3] describes the genera- 
tion of grids of model SEDs and describes the fitting 
procedure by means of which the observational data are 
compared to these SED grids to identify the best-fitting 
model. The detailed derivation of the fitting formalism, 
including the formulae for the case when some of the pho- 
tometric observations yield only non-detections, is given 
in the Appendix. Section [5] gives several examples of 
the application of the SEDf it software and S}6] discusses 
several issues associated with SED fitting in general. 

2. OVERVIEW AND IMPLEMENTATION 

Given a distant galaxy, the task at hand is to compare 
a set of photometric observations made through (usually) 
broadband filters with sets of flux densities predicted for 
a grid of spectral models, in order to identify that model 
which best matches the data and to establish the degree 
of confidence in that model. The task naturally sepa- 
rates into two distinct components: (1) the generation of 
model flux densities, and (2) the quantitative comparison 
of these model flux densities to the observational data. 
SEDfit deals with these two components by means of 
two separate software programs: make_sed and f it_sed 
(Fig.P). 

The program make_sed carries out the tasks associated 
with the first of these two major components, namely the 
generation of a grid of model magnitudes (or spectra) . It 
starts by taking as input a set of rest- frame model spectra 
(from, e.g., the GALEXEV library of Bruzual & Chariot 
2003) and performs the following operations: attenua- 
tion by interstellar dust, cosmological dimming and red- 
shifting to the observed frame, and redshift-dependent 
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attenuation by intervening intergalactic gas clouds (see 
top half of Fig. [1]). The spectra thus adjusted can then 
be saved, or can be further integrated through desired 
system transmission curves to produce grids of model 
magnitudes. The details of this process are described 
in SjHl The program make_sed loops over the full speci- 
fied ranges for all the user-specified parameters (redshift, 
reddcnning, etc.), as well as parameters internal to the 
input spectral templates (such as stellar population ages) 
to produce grids of models as a function of these param- 
eters. 

Next, the program fit_sed takes the grid of model 
magnitudes generated by make_sed and compares it to 
the observed photometry to find the model that best 
matches the data. The comparison is done in flux space 
(after converting from magnitudes to flux densities) by 
means of a maximum likelihood test ( "x 2 fitting" ) mod- 
ified to correctly account for those photometric obser- 
vations that have only yielded upper limits. The de- 
tails of the mathematical formalism that underlies the 
fitting procedure used by f it_sed are given in <2J In ad- 
dition to identifying the best-fitting model and provid- 
ing the corresponding parameters and their uncertain- 
ties, f it_sed can optionally also produce \ 2 maps of the 
parameter space, Monte Carlo simulations of the uncer- 
tainties, and/or a means to produce the model spectra 
that correspond to the best-fitting model. The program 
can loop over multiple objects in a catalog, fitting each 
one with the full grid of models or a subset thereof. 

The present implementation of the software is in Perl 
and the Perl Data Language (Glazebrook & Economou 
1997). While Perl itself is a scripting rather than a com- 
piled language, the Perl Data Language module is very 
efficient at dealing with large arrays because it uses pre- 
compiled libraries. The computationally-intensive work 
in SEDf it is done using this Perl Data Language func- 
tionality and is therefore quite fast. In particular, array 
operations are very efficient in Perl Data Language, mak- 
ing it possible for SEDfit to perform much of the SED 
generation and fitting very rapidly. 

Both make_sed and fit_sed are operated from the 
Unix command line and are controlled by means of 
editable parameter files and command-line overrides of 
those files. This approach and its syntax is modeled 
on those adopted by the SExtractor photometry package 
(Bertin & Arnouts 1996) and will be intuitively familiar 
to those experienced with SExtractor. This approach al- 
lows for fast experimenting with parameter settings and 
enables efficient scripting. 

3. GENERATION OF MODEL FLUXES AND MAGNITUDES 

The process of generating model SEDs usually begins 
with unattenuated rest-frame spectra. These input spec- 
tra can be either empirical (e.g., Coleman, Wu, & Weed- 
man 1980) or synthetic (e.g., Bruzual & Chariot 1993, 
2003) and can represent galaxies, AGN, individual Galac- 
tic stars, or other objects. SEDfit does not provide a 
built-in ability to mix spectra to, e.g., create more com- 
plex star formation histories than those already provided 
in the spectral libraries or to add emission lines to spec- 
tra, such as the Bruzual & Chariot (2003) library, that 
lack them (but see, e.g., Yabe et al. 2009). It is left up to 
the user to combine or modify spectra if so desired and 
then save them in the format readable by SEDfit. This 



approach allows maximum flexibility in what models can 
be fit and SEDf it's only requirements are that the in- 
put spectra follow a standardized format and are given 
in units of power per frequency interval as a function of 
wavelength in angstroms, L V {X). This system of units 
reflects the tendency of high-z optical/IR astronomy to 
use AB magnitudes (which are defined as logarithms of 
power per frequency intervaQ) but to work as a function 
of wavelength rather than frequency. Then, starting from 
L„(A), several steps are needed to arrive at the spectral 
energy densities in the observer's frame, f v (X). These 
steps are performed by the SEDfit program make_sed as 
described below. 

3.1. Attenuation by interstellar dust 

First, each input spectrum L v (\) is modified to simu- 
late the effect of attenuation due to interstellar dust in 
the model galaxy. This attenuation can be expressed as 

LUA) = MA)-10-°- 4S(B - v)fc(A) , (1) 

where L' and L are the attenuated and unattenuated 
spectral energy distributions, E(B—V) is the color excess 
whose value controls the overall amount of attenuation, 
and k(X) is the dust law which defines the profile of atten- 
uation as a function of wavelength. At present make_sed 
supports the following dust laws: the Calzetti (1997) and 
Calzetti et al. (2000) starburst dust laws, the Fitzpatrick 
(1986) Milky Way, LMC, and 30 Doradus dust laws, and 
the Prevot et al. (1984) SMC law. The user specifies one 
of the dust laws and one or more E(B — V) value(s) to 
be included in the calculation of the model grid. 

3.2. Attenuation by intergalactic gas 

Next, the spectral energy distribution is attenuated 
further to account for the stochastic effects of absorption 
due to intergalactic gas lying between the object being 
simulated and the observer. At present two intergalac- 
tic attenuation prescriptions are included in make_sed: 
Madau (1995) and Inoue et al. (2005). The attenuation 
operation can be expressed as 

L'X\) = L' v {\).e- r '"< x >'\ (2) 

where L" is the SED attenuated by intergalactic gas, 
L' is the dust-attenuated SED of eq. [TJ and r e t / is the 
effective optical depth as specified by Madau (1995) or 
Inoue et al. (2005). This optical depth, r e //, depends on 
the redshift of the source because the absorbers' column 
density increases with distance and the absorption wave- 
lengths are subject to the cosmological redshifts of the 
absorbers. 

The user can chose one of the two intergalactic atten- 
uation prescriptions and also specify one or more atten- 
uation strengths. The default absorption strength rep- 
resents a typical line of sight through the Universe, but 
sightlines that are less or more crowded can also be sim- 
ulated. Additionally, since the attenuation is redshift- 
dependent, one or more redshifts — representing the 
rcdshift(s) of the object(s) being simulated — need to 
be specified. 

1 Recall that m A B = -2.5 log [f v / (erg s^cm^Hz- 1 ] - 48.60 
(Oke 1974), or m AB = -2.5 log(/„/nJy) + 31.4. 
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Fig. 1. — Flowchart of SEDfit operations. The upper part of the diagram illustrates the operations performed by make_sed, the program 
that generates model spectral energy distributions. The lower part shows the operations performed by f it_sed, the program that performs 
the comparisons between models and data. 
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3.3. Cosmological redshifting and dimming 

The next step is to apply the effects of cosmological 
distance to the spectrum and thus transform it into the 
observed frame flux density distribution. This operation 
can be expressed as 



L'UX/jl + z)) 
AnDl/il + zY 



(3) 



where Dl is the luminosity distance. The A/(l + z) ex- 
pression in the numerator indicates that the spectrum's 
wavelength scale has been redshifted from the original, 
rest-frame wavelength to one at the desired redshift. The 
user can chose a desired cosmology by specifying the val- 
ues for f2jvf , £2a , and Hq , while the same single or multiple 
redshift values are used as for the intergalactic attenua- 
tion of eq. [2] 

3.4. Generating model magnitudes 

As the result of the above steps, make_sed contais a 
grid of redshifted and attenuated spectral density dis- 
tributions. These SEDs can be saved to a file as a 
grid of spectra, or can be further processed to generate 
broadband SEDs by averaging (sometimes called "con- 
volving" in the field) the spectra through system trans- 
mission curves, which can be provided by the user and 
can/should include the filter transmission profile, detec- 
tor efficiency, and other wavelength-dependent through- 
put characteristics of the telescope+instrument system. 
For photon-counting systems, such as CCDs and near-IR 
arrays, the procedure can be expressed as 



fn 



JM\)\T t (X)d\ 
J XTi(X)dX ' 



(4) 



where Tj is the system transmission curve for the i 
filter and f m i is the model's observer-frame flux density 
through that filter (the subscript m refers to "model" 
as distinct from the observed "data" flux densities to be 
encountered in @ . The A factor in both numerator and 
denominator accounts for the fact that the number of 
photons per unit flux varies with wavelength. The user 
can specify one or more transmission curves to be used; a 
large number of common transmission curves is already 
supplied with make_sed and more can be generated by 
the user as simple ASCII files. If the input spectrum L„ 
is in units of erg/s/Hz then the output of make_sed is in 
AB magnitudes, 



m A B. 



-2.5 log /„ 



.60, 



(5) 



(Oke 1974). 



3.5. Outputs 

The main end product of make_sed is a file that con- 
tains one galaxy model per line, with each model con- 
sisting of a number of m,i magnitudes along with associ- 
ated parameters such as redshift, E(B — V) value, galaxy 
stellar mass, etc. The number of models in the grid is 
Nmodeis = N E{B _ V) x N z x N spectra , where N E(B _ V) is 
the number of color excess steps sampled, N z is the num- 
ber of redshift steps, and N spec t ra is the number of input 
spectra (the F v ) that were used. These input spectra 
carry with them information such as star formation his- 
tory, metallicity, etc., and the user can provide as few or 



as many input spectra as are desired. A typical make_sed 
output file may include 221 spectra (the standard num- 
ber provided for a given metallicity and star formation 
history by the GISSEL package of Bruzual & Chariot 
2003), several hundred redshift steps, and several tens of 
E{B - V) steps. 

This grid of models produced by make_sed can be used 
for exploring the expected magnitudes and colors of high- 
redshift objects, but more pertinently here it serves as 
one of the two ingredients (the other being the obser- 
vational data) in determining the redshifts or physical 
properties of observed high-z galaxies, as described in 

m 

In addition to the model magnitudes grid, make_sed 
can produce a file that contains the full redshifted and 
attenuated spectra before the integration through system 
transmission curves (i.e., the end product of £|3.3[) . Such 
spectra are not used directly in SED fitting (Q but are 
useful for, e.g., visualizing the SEDs of distant galaxies, 
as is illustrated in the left panels of Figs. [2] and [3] (©. 

4. SED FITTING 

SEDf it uses a maximum- likelihood approach to deter- 
mine the best-fitting model and to establish confidence 
regions in parameter space. These tasks are performed 
by the program fit_sed. As its inputs fit_sed takes 
the magnitudes model grid produced by make_sed (Sj3j 
and a catalog file containing the observed magnitudes or 
flux densities. The format of the catalog file is flexible 
as long as it contains one object per line. Operationally, 
fit_sed determines the best- fitting model for a given 
object by comparing the object's flux densities with the 
model flux densities which are obtained from the model 
grid file. The model that is the most likely, in the sense 
that it gives the lowest % 2 value, is deemed to be the 
best-fitting model. 

4.1. Determining the goodness- of- fit for a model 

When fitting an object, f it_sed properly accounts for 
both detected and observed but undetected fluxes. When 
all the photometric data are detected flux measurements, 
the f it_sed compares a model with the observed data 
and determines the quality of the fit by means of the 
usual x 2 statistic, 



2 ( fd,i ~~ s .fn 



(6) 



where f m ^ is the model flux density through the i th band- 
pass, fd,i is the observed flux density that same bandpass, 
and o~i is the uncertainty in the observed flux density. 
Note that this uncertainty could include both the claimed 
observational photometric error and an additional term 
that captures an estimate of systematic uncertainty. 
SEDf it does in fact allow the user to specify a ai >sys for 

each bandpass, in which case <7, = {o~i. p hot + o'i.sys) 1 ^ 2 . 

The quantity s in Eq. |6] is the flux scaling between 
the model and the data which is determined analytically 
(Sawicki 2002) as 



fd,ifrr 



E 



f 2 

J n 



(7) 
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The derivations of eq. [5] and [7J are reviewed in the Ap- 
pendix. 

Equation [6] yields the correct result when the object 
is detected in all the observed bandpasses but becomes 
ambiguous when one or more of the observed band- 
passes yield such a non-detection. However, even non- 
detections can be useful or even critical (identification 
of Lyman break galaxies provides an obvious example in 
this respect) in constraining models or ruling out por- 
tions of parameter space. Consequently, rather than ig- 
noring the information offered by such non-detections, 
or adopting the ad-hoc approach of assigning arbitrary 
fluxes to the non-detections, the current implementation 
of SEDfit properly accounts for such data by using a 
modified version of the \ 2 method. The derivation of 
this modified \ 2 statistic is given in the Appendix and 
yields 



fd.i $frt 



2£ln 



flin 



exp 



/-*/, 



df. (8) 



The first sum on the right-hand side in eq. [8] refers to 
detected bandpasses, indexed with i, while the second 
sum accounts for those bandpasses, indexed with j, that 
yielded upper limits; here fu m ,j is the 1-a flux detection 
treshold in the j th such bandpass, which f it_sed calcu- 
lates from the n-a values (where n=3, 5, etc.) given in 
the data catalog. Note that when all bandpasses yield 
detections, the second term in eq. [8] drops out and eq. |8] 
reduces to the familiar eq. El as expected. Whereas the 
flux scaling s could be determined analytically for the 
case when all bandpasses yielded detections (eq. [7J , no 
such exact analytic expression exists for the case when 
non-detections are present. For this reason, the integral 
in eq. [8] needs to be evaluated numerically, making the 
fitting of objects with non-detections significantly more 
computationally expensive than that of objects that are 
detected in all bandpasses. Nevertheless, non-detections 
can be critical in ruling out parts of parameter space (see 
§ 15.21 for an example) and thus the additional computa- 
tional cost is often worthwhilfl 

4.2. Identifying the best-fit model 

The program f it_sed searches the model grid to find 
the model that minimizes x 2 (as defined by eq. |8|) . Two 
search algorithms are implemented in the standard re- 
lease of f it_sed and the user is free to choose one or the 
other. 

The simpler and most robust, "brute- force" fitting ap- 
proach computes \ 2 f° r every model in the grid and 
then unambiguously identifies the one with the lowest 
X 2 value. This approach is robust but — because of the 
need to calculate % 2 for every model in the grid — not 
particularly efficient. However, in many applications it is 

2 Were all flux measurements reported, even those below a nom- 
inal detection threshold, it would be possible to use eq. [7] in all 
situations. However, it is usually the case in astronomy that only 
limits are reported when flux measurements are below an assigned 
threshold, and the approach taken by eq. \E\ provides a way to in- 
clude that information in the fit. 



sufficiently fast while its robustness against effects such 
as secondary minima as well as its ability to produce 
complete \ 2 maps of the model space are its advantages. 

The second available option is to perform a downhill 
search of the model parameter space starting from a ran- 
dom position in the model grid. Because the downhill 
search can become lodged in local minima in parameter 
space, fit_sed can repeat the downhill search several 
times per object, starting with a new random position 
each time. The user specifies the number of searches 
per object; experiments show that ~10 repeats per ob- 
ject are typically sufficient to ensure that the absolute 
minimum is found. Because of its speed, this "downhill 
search with random repeats" option is preferable when 
extremely large sets of objects are to be fitted, or - 
because of the computationally-expensive nature of eval- 
uating non-detections (eq.[8]) — when many objects with 
upper limits are present. 

Additional search algorithms (such as simulated an- 
nealing) have been implemented within fit_sed but it 
was found that the two standard approaches (the brute- 
force calculation for every model in the grid and the 
downhill search with random repeats) are adequate in 
the vast majority of situations; consequently at present 
only these two algorithms are included in the standard 
distribution of SEDfit. 

The model with the lowest \ 2 value is by definition 
the most likely to have produced the observed data and 
is adopted as the "best-fitting" model. The values of the 
parameters associated with this model (such as redshift, 
age, spectral type, etc.) are then adopted as the best- 
fitting parameter values. 

4.3. Determining the uncertainties 

Uncertainties on the values of parameters in the best- 
fitting model can be established by fit_sed in one of 
several ways. 

The first way is to make use of the Monte Carlo simu- 
lations feature of f it_sed. Here, f it_sed refits each ob- 
ject a large, user-specified number of times (e.g., 200 or 
1000), but with each realization the photometry in each 
band is randomly perturbed from the actual value by 
drawing from a Gaussian distribution with standard de- 
viation dictated by the object's photometric uncertainty 
in that band. Note that for operational convenience the 
perturbation is performed on the model photometry; this 
is mathematically equivalent to perturbing the observed 
data and is consistent with the fundamental principle 
that underlies the \ 2 formalism. The distribution in pa- 
rameter space of the best fits to the perturbed photom- 
etry realizations then defines the uncertainty ranges on 
the parameters of the true best- fit model. SEDfit reports 
the uncertainties on each parameter in its tabulated out- 
put, and the individual Monte Carlo results can also be 
saved to a file for establishing confidence contours ( "error 
contours" ) in multi-parameter space. The Monte Carlo- 
based uncertainties work when all the photometric bands 
are detected as well as in the case when some bands are 
undetected. 

An alternative, faster way, that can be used when the 
full grid of models has been fit (i.e., the "brute-force" 
fit, ij4.2p is to use the \ 2 map produced during the fit- 
ting process. The offset between each model's \ 2 value 
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0.5 1 1.5 2 0.2 0.4 0.6 1 2 3 4 5 

A obs ( M m) E(B-V) SFR(M /yr) 

Fig. 2. — An example of SED fitting of a Lyman Break Galaxy in the Hubble Deep Field. This faint galaxy (i?,=27.05) has no known 
spectroscopic redshift (redshift is a free parameter of the fit) and in this example has been fit with a suite of solar metallicity, constant star 
formation rate models from the Bruzual & Chariot (2003) library. The left panel shows the photometric data points and the best-fit model 
spectrum (shown as a solid line), with the best-fit model's parameters listed in the figure. The middle and right panels show the location of 
the best-fit model in parameter space along with an indication of the associated uncertainties. The middle panel shows confidence regions 
determined from x 2 maps: the darker shading indicates higher \ 2 values and the contours correspond to 1 and 2 a confidence. The right 
panel shows the best-fit model along with 1000 Monte Carlo re-fits of the data; the distribution of these 1000 realizations gives an indication 
of the probability distribution of the model parameters; although this is not done here, formal confidence regions could be generated from 
these Monte Carlo results by converting them into a two-dimensional histogram and plotting density thresholds. 



from the x 2 of the best- fitting model, Ax 2 , is a measure 
of probability of the corresponding model. On the basis 
of this offset the one-dimensional uncertainties on each 
parameter as well as full confidence contours in multi- 
parameter space can be established. The A% 2 level that 
corresponds to a desired confidence level for the number 
of free parameters available can be taken from theory 
(e.g., Press et al. 2007) or can be established empiri- 
cally for a given dataset using a representative sample of 
Monte Carlo simulations. 

4.4. Outputs 

The program f it_sed produces one standard and sev- 
eral optional output files. The standard output file 
lists, for each object in the input catalog, the best- 
fitting model and its associated parameters and their 
one-dimensional uncertainties. The user can specify 
which best-fit model parameters to report in the stan- 
dard output file and also which parameters from the in- 
put data catalog to propagate into that file. The three 
optional outputs arc: a file that contains the fit results for 
all the Monte Carlo realizations; a x 2 map for each object 
in the input file; and a file that contains the shell-style 
command-line input which, when supplied to make_sed, 
will result in the production of the redshifted and dust- 
attenuated spectra that correspond to the best-fit mod- 
els. All these optional outputs can be used to illustrate 
the quality of the SED fit; some examples of the uses of 
such outputs are given in the figures in SJS] 

5. EXAMPLES OF APPLICATIONS 

This section illustrates the use of SEDf it by means 
of three examples: (1) the basic SED fitting of a high- 
z galaxy, (2) SED fitting that makes use of photomet- 
ric non-detections, (3) two-dimensional SED mapping of 
high-z galaxies, and (4) simulations of photometric red- 
shift performance. Of course, SEDf it is not restricted to 
just these four specific tasks and many additional exam- 
ples of SEDf it applications may be found in the papers 
listed in fJTJ 



5.1. Simple SED fitting of a distant galaxy 

The determination of quantities such as stellar mass, 
extinction, age, etc. constitutes SEDf it's main raison 
d'etre. As an illustration, an example of the results of the 
fitting of a high-z galaxy is given here. Figure[5]shows the 
photometric properties of a faint color-selected Lyman 
Break Galaxy in the Hubble Deep Field and illustrates 
various aspects of SEDf it's output, including the best- 
fitting spectral model (left panel) and confidence regions 
in parameter space (middle and right panels) determined 
in two different ways. 

Although, for simplicity, Fig. [2] presents the results for 
a single galaxy, SEDf it is designed to routinely process 
not just individual objects but catalogs of galaxies. Fur- 
ther examples of such work can be found in Sawicki & 
Yee (1998), Yabe et al. (2009), and Sawicki (2012). 

5.2. SED fitting including non- detections 

SEDf it's ability to include non-detections in the fitting 
process is illustrated in Fig. [3J In this example the issue 
at stake is constraining the age of a stellar population of 
a z=1.5 object. The left panel of Fig. [3] shows both ob- 
servations that yielded detections (BVI) and those that 
resulted only in upper limits (JH). 

When only the detected points are used in the fit, the 
data altogether fail to constrain the age of the object: 
the full range of ages between very young (< 10 7 yr) and 
that of the age of the Universe at the object's redshift 
is allowed, as is shown by the gray contour in the right 
panel of Fig.[3j In contrast, when the non-detections are 
included in the fit (using the full eq. [5]), a more useful 
constraint on the age is possible: at the 68% confidence 
level, the age is younger than 5 x 10 s yr (black contour). 

The formal maximum- likelihood fit (eq. [5J gives a sta- 
tistically correct way to incorporate the information con- 
tained in the non-detections in a way that a "%-by-eye" 
approach of inspecting the models in the left panel of 
Fig. [3J does not allow. The fact that some permissible 
models (black curves) in the left panel of the figure ap- 
pear higher than the ler JH limits is simply because such 
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Fig. 3. — An example of SED fitting that uses non-detections to constrain the fit. The object is a galaxy with redshift fixed at z = 1.5, 
photometric detections in the BVI bands and upper limits in J and H (indicated in the left panel with downward-pointing arrows at the 
1<t flux density limits) fit using Bruzual & Chariot (2003) constant SFR models. A sampling of models that give acceptable fits to the 
detected data (i.e., BVI only) is shown with light-colored curves in the left panel. Models that give acceptable fits once non-detections are 
included in the fit are shown with black curves. The right panel shows regions of &ge-E(B — V) parameter space that are permitted by the 
fits: the light contour encloses models allowed when only detections are used in the fit, while the black contour shows the constraints when 
non-detections are included. The contours show 68% confidence regions in both cases. The inclusion of non-detections clearly improves the 
constraints on the parameter space, here enabling the placement of meaningful constraints on the age of the stellar population. 



models are indeed permitted by the data. The most 
tempting "x-by-eye" analysis here would have missed 
these older but permissible models to give very low but 
unrealistic constraints on the age. In contrast, the proper 
maximum-likelihood treatment that uses Eq. [S] gives 
constraints that are both more accurate, automatic, and 
reproducible. 

5.3. Spatially resolved SED fitting 

SEDfit can also be used to perform spatially-resolved, 
pixel-by-pixcl SED fitting when the target galaxies are 
sufficiently resolved in the images. Operationally, the 
procedure is a simple modification of the procedure used 
to fit a catalog of individual objects f ^5.1[) . In this ap- 
proach, each pixel of a galaxy is fit as if it were an indi- 
vidual object and then the resulting lists of best-fit values 
are reassembled into parameter maps (e.g., Abraham et 
al. 1999). 

Figure |4] illustrates the results of this procedure ap- 
plied to two Lyman Break Galaxies in the Hubble Deep 
Field. Here, HST l%oo-B45oVk)6-f8i4</iio-ffi60 images of 
spectroscopically-confirmed z^2.3 galaxies have been 
smoothed to a common PSF to guard against artificial 
color gradients, their individual pixels fit using SEDfit, 
and the results reassembled into spatial maps of SFR, 
stellar mass, color excess E(B — V), and age (Palmer 
2010). The resulting maps reveal interesting spatial off- 
sets between regions of ongoing star formation and more 
quiescent regions of older stars. 

5.4. Simulation of photometric redshifts 

SEDfit can be used to estimate photometric redshifts 
to distant galaxies (e.g., Sawicki et al. 1997; Sawicki & 
Mallen-Ornelas 2003; Sorba & Sawicki 2010) and also 
to perform simulations of photometric redshift perfor- 
mance (e.g. Sorba & Sawicki 2011). Figure [5] shows an 



example of such simulations: the specific case here shows 
the expected fidelity of the photometric redshift perfor- 
mance of the proposed Euclid dark energy mission when 
augmented by the addition of [/-band photometry to the 
planned Euclid filter set. In this example, photometry for 
several hundred thousand simulated galaxies was gener- 
ated using make_sed, drawn from a range of spectral tem- 
plates and using a realistic redshift distribution. These 
galaxies' photometry was then perturbed according to 
the expected photometric uncertainties. Finally, the "ob- 
served" galaxies were fit using fit_sed with the same 
spectral templates as used in their generation. Figure [5] 
shows the "observed" photometric redshift vs. the input 
"spectroscopic" redshifts; because the same spectral tem- 
plate set was used to generate the artificial galaxies as 
to SED-fit them, the scatter in the "observed" redshifts 
from the diagonal line results purely from photometric 
uncertainties in the data. Analysis of this type can be 
useful when designing photometric redshift surveys. 

6. DISCUSSION 

Given observational photometric data and a set of 
model spectra, SED fitting can produce useful insight 
into the properties of distant objects. In this regard, the 
most commonly-sought parameters are photometric red- 
shifts, galaxy stellar masses, amount of interstellar ex- 
tinction, dust-corrected star formation rates, and ages of 
stellar populations. The technique has now been widely 
used and forms a standard part of the toolbox of extra- 
galactic astronomy. 

Despite the apparent simplicity of SED fitting and the 
clear usefulness of the physical properties that it purports 
to derive from pure photometric data, it is important to 
use the technique with caution since there are a num- 
ber of issues that can bias the results. Some of these 
limitations are discussed in the following paragraphs. 
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Fig. 4. — Best-fit parameter maps for two Hubble Deep Field 
z~2.3 Lyman Break Galaxies whose individual pixels were fit using 
SEDf it. Each panel is 30 pixels, or 1.2", on the side, which corre- 
sponds to ~10 kpc at these redshifts. These maps suggest that in 
both galaxies unobscured older stars are offset from younger, more 
dust-obscured regions of recent star formation (Palmer 2010). 



First, no matter how good the spectral models, the 
observational data need to be of the right quality. This 
means not only a sufficiently large and well-chosen set 
of photometric bandpasses, but also accurate photome- 
try, since otherwise systematic errors in, e.g., photomet- 
ric calibration, will propagate into systematic errors in 
the estimated parameters such as stellar masses and star 
formation rates. In this respect it is crucial to realize 
that accurate colors (flux density ratios) are far more 
important than the absolute calibration of the overall 
spectrum: for example, a 10% offset in overall, across- 
all-bands calibration will result in only a 10% offset in 
stellar mass and star formation rate, but a 10% offset 
in relative photometry between bands can result in far 
worse and even catastrophic inaccuracies in the fitted 
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Fig. 5. — Photometric redshift simulations for the Euclid dark 
energy mission performed using SEDf it. The intensity of color in 
the figure corresponds to the number of galaxies in each z pflot - 
Zspec bin, with brighter colors corresponding to more objects. The 
numbers in the lower right of the figure indicate (top to bottom) 
number of galaxies per square arcminute, the median redshift, the 
photometric redshift bias, and the fraction of catastrophic redshift 
failures, de fined here as the fraction of objects outside the dashed 
lines (see £15.41 and Sorba &; Sawicki 2011 for additional details). 
Analysis of simulations like these provides the ability to fine-tune 
the design of photometric redshift surveys. 

quantities. Careful matching of photometric apertures 
across bandpasses to ensure accurate flux density ratios 
is thus essential, including procedures to account for PSF 
differences between bands. 

A second set of limitations stems from the inability of 
SED fitting to constrain certain key properties. A well- 
known limitation (e.g., Sawicki & Yee 1998; Papovich et 
al. 2001) is the inability to constrain star formation his- 
tories from the typical broadband photometry available 
for distant objects. The degeneracy in star formation his- 
tories propagates into degeneracies in other fit parame- 
ters, most notably the star formation rates. There is 
no good way around this issue given limited photomet- 
ric data and thus one simply has to be aware of this 
limitation when interpreting fitting results. In this re- 
spect it is perhaps wise to perform relative — rather 
than absolute — measurements, comparing results across 
similarly-selected samples of galaxies fit with a single set 
of models rather than trusting the absolute results of the 
best-fit parameters. 

A further limitation is the SED fitting technique's fo- 
cus on the most luminous stellar population present in 
the object being fitted. This limitation most commonly 
manifests itself in the inability to properly account for 
an older stellar population when a young population is 
present. The limitation stems from the fact that old stel- 
lar populations are faint in rest-frame UV compared to 
young populations that are dominated by luminous high- 
mass stars. Old stellar populations, if they exist, will 
thus often be missed in the glare of even a sprinkling of 
bright young stars. This effect can be critical since such 
old, undetected underlying stellar populations can con- 
tain significant masses of stars, potentially several times 
the mass of the young detected populations. 
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An additional limitation stems from the input spec- 
tra used to generate the model templates: the results of 
the SED fitting will be only as good as these spectra. 
While a number of excellent stellar population models is 
now available, significant differences between these mod- 
els remain and suggest that our understanding of stellar 
populations is still not perfect. The resulting estimates 
of galaxy physical parameters will thus carry associated 
systematic errors. In light of such differences it is per- 
haps wise to again rely on relative comparisons between 
fit galaxy properties (masses, star formation rates, etc.) 
rather than on the absolute values of these properties. 

A related point is that SED fitting will produce a best- 
fitting model regardless of whether that model is or is not 
a good match to the observations. It is therefore wise to 
inspect the goodness-of-fit parameter, % 2 , to ensure that 
the best-fitting model is indeed a well-fitting model (the 
two are not synonymous!). Further, irrespective of how 
good a x 2 is produced, even a well-fitting best-fit model 
may not be a true reflection of reality since degenera- 
cies in fitting space exist and the SED-fitting procedure 
may preferentially choose a slightly better-fitting, but ul- 
timately incorrect model over a worse-fitting but correct 
one. Consideration of such degeneracies should be an 
important aspect of interpreting SED-fitting results. 

Despite systemic limitations such as those discussed 
above, SED fitting has proven to be a useful and pop- 
ular technique in extragalactic research. It provides us 
with a way to estimate key quantities of interest, such 
as photometric redshifts, galaxy stellar masses, amount 
of extinction due to dust, dust-corrected star formation 



rates, and ages of stellar populations; and it lets us do so 
for a large number of objects from a single set of common 
observations. The technique has served us in this role for 
over a decade now and will likely continue to be a key 
tool in our toolbox well into the coming era of ever more 
powerful observational facilities and ever larger datasets. 
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ment the Inoue et al. (2005) intergalactic opacity func- 
tion in the present version of the software; to Michael 
Palmer, Bobby Sorba, Kiyo Yabe, and Suraphong Yuma 
for requests that resulted in extending the capabilities of 
this software; to Bobby Sorba for allowing me to show 
the results of our Euclid photo-z simulations in Fig. [5] 
and to Michael Palmer for letting me display results 
from his undergraduate honours thesis in Fig. 2J I also 
thank Liz Arcila-Osejo, Anneya Golob, Michael Palmer, 
Taro Sato, Barbara Sawicki, Jerzy Sawicki, Bobby Sorba, 
Kiyo Yabe, and Suraphong Yuma and the anonymous 
referee for comments that improved the quality of this 
manuscript. Finally, this software has evolved over a 
number of years, and thus funding from several agen- 
cies has contributed to its development: I gratefully ac- 
knowledge support from the Natural Sciences and En- 
gineering Research Council (NSERC) of Canada, the 
Canadian Space Agency (CSA), NASA, an International 
Long-Term Visitorship from the Japan Society for the 
Promotion of Science (JSPS), and the Atlantic Compu- 
tational Excellence Network (ACEnet). 



APPENDIX 



THE MAXIMUM LIKELIHOOD FORMALISM FOR SED FITTING WITH UPPER LIMITS 

SEDf it's model-fitting program, f it_sed uses a maximum-likelihood approach to identify the best-fitting model and 
determine confidence regions in parameter space. This appendix first reviews the derivation of the usual % 2 formula 
as applied to the particular case of SED fitting and then uses a similar approach to derive a variant of \ 2 that is 
appropriate when some of the data are non-detections. 

The standard way to approach the development of a maximum-likelihood fitting method is to assume that a given 
model underlies the observed data and the data are drawn from that model but are perturbed by Gaussian errors (see, 
e.g., Press et al. 2007). In other words, the question one asks is "what is the probability that the observed data could 
have occurred if the data were drawn from the model that's being considered?" 



The case when the object is detected in all bands 

If the object has been detected in all photometric bands then the approach outlined above produces the well-known 
X 2 formalism, as is shown below (see, e.g., §15.1 of Press et al. 2007, for a general x 2 derivation). For a single band i, 
the probability that a set of observations has been produced from a given model is 



Pi oc exp 



1 ( fd,i - Sf„ 

2 



A/, 



(Al) 



where fd,i is the observed (data) flux density in the ith band, o~i is its standard deviation (i.e., uncertainty), and s/ m> j 
is the model flux density in the same band. The quantities involved are illustrated in Fig. Ela). Note that we have 
separated the model flux density into two components: the quantity s is the normalization of the model, which, for 
a given model, is the same across all the bands i, while f m ^ contains the information about the spectral shape of 
the model. This separation allows us to considerably speed up the numerical handling of the calculation, as is shown 



below (eq. IA5|) . In practice, having f m ^ distinct from the normalization s is also convenient because it is f m ^ that is 
produced from the spectral models (eq. [4]) , while s contains the information about the normalization adjustment that 
best brings such a model into match with the observed data and gives the multiplicative scaling that links the model's 
scaleable parameters — such as luminosity, stellar mass, and star formation rate — with those in the galaxy under 
consideration. 

The overall probability of the data in all the observed bandpasses being drawn from the model in question is then 
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(a) detection 




f d,i sf m,i 
(b) non— detection 




Fig. 6. — Probability that an observation is drawn from a given model. The top panel is for the case of a detection, which leads to 
the well-known standard x 2 goodness-of-fit estimator. The bottom panel is the equivalent diagram for the case of a non-detection. In 
both cases the model predicts the flux density value sf m and the Gaussian curves show the probability distribution under the influence 
of Gaussian perturbation of the flux. The shaded regions (in practice infinitessimally thin in the case of a detection) correspond to the 
probability that the random Gaussian uncertainties have perturbed the true flux to result in the observed detection at j (top panel) or 
non-detection below the detection threshold fu m j (bottom panel). 



proportional to the product of the individual probabilities, 

PocJJP,. (A2) 



Next, taking the logarithm of eq. IA2I gives 



P oc In J] P t = - \ g ( h ' 1 a S fm ^ ) 2 + g ^ A/. (A3) 



The second sum in eq. IA3I is a constant and thus maximizing the probability is just equivalent to minimizing the first 
sum, 

X^Ef ^V^'O 2 - (A4) 

The derivation thus far has followed the standard derivation of the % 2 formula as a maximum likelihood estimator 
with the assumption of normally distributed errors (e.g., Press et al. 2007). 

The normalization factor s of the best-fitting model can be derived analytically and doing so can considerably speed 
up the numerical calculation of % 2 . This analytic solution can be achieved by minimizing the x 2 of eq. IA4I by taking 
the partial derivative dx 2 /ds and setting it to zero (Sawicki 2002). Doing so yields the expression 

Efd,ifm,i I fm,i / a r\ 

/2^^ 2 "' (A5) 

i % i 1 

which exactly gives the appropriate scaling s of the model. Identifying the model with the smallest \ 2 from among the 
model suite gives the most likely (i.e., best- fitting) model in the case when the object is detected in all the observed 
bands. 

The case of upper limits 

For the case where the object is not detected in one o r more of the observed bands we can proceed as above but 
with replacing detections with upper limits in eq. IAH - |A"4l In the case of a non-detection in the band j, the probability 
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that the observation in that band is drawn from a given model is 



exp 



1 ( f - Sfm, j 

2 



df, 



(A6) 



where fum,j is the upper limit (i.e., sensitivity) of the observation in the jth band. In analogy to eq. lA2[ the probability 
of the data (detections and non-detections) across all the observed bandpasses being drawn from the given model is 
then the product of the individual probabilities, 



(A7) 



The subscripts i indicate those bands which have yielded detections, while the subscripts j are for the bands with 
non-detections. The equivalent of eq. IA3l then is 



P ex 



1 



E 



fd.i sfrr 



^lnA/ + ^ln/ "exp 



I ( f - S f m, 3 
2 



df, 



(A8) 



and, in analogy with eq. IA41 maximizing the likelihood that the observed dataset (detections and non-detections) is 
drawn from the given model is equivalent to minimizing the quantity 



E 



fd.i sfn 



2£ln 



/fir, 



exp 



1 ff-sf; 

2 



m.,0 



df. 



(A9) 



No te th at if all the observed bands yielded detections then the last term drops out and eq. IA9I reduces to the form of 
eq. IA4| as one would expect. 

For computational convenience, the integral in eq. IA9I can be recast in terms of the error function, erf (a;) = 



(2/0F) fie'* dt, so that 



x 2 = 



E 



fd,i sfri 



- 2 !>{v! ff ' 



1 + erf 



flim,j ^fm,j 



(A10) 



Finding the most likely model requires evaluating eq. lAlOl to identify Xmini the smallest value of x 2 - I n the case when 
all the bands yield detections finding Xmin can accelerated by optimizing the model scaling factor s using eq. IA5I 
In the present case of non-detections there is no simpl e equi valent of eq. IA5I One approach is then to find the s that's 
optimal for a given model by numerically exploring eq. IA101 Alternatively, in analogy with eq. lA51 setting d\ 2 /ds = 
gives the condition for the optimal s for a given model: 



E 



Mi - sf m A ff m ,i\ ft /mjexp {" ^ flim,j ~ sfm <^ a ^ 2 } 



Jm,i \ / * \ ^ 



"Kj G i { 1 + eri [{flitnj ~ sfm,j)/V2crj] } 



(All) 



Finding the root of eq. IA11I then gives the optimal s for the model. This root can be obtained numerically using a 
variety of root-finding methods. The present implementation of SEDfit uses simply numerically searches for a model's 
optimal s in eq. IA10I 
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